home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / newmat.exe / EXAMPLE.CXX < prev    next >
C/C++ Source or Header  |  1991-07-30  |  7KB  |  178 lines
  1. //$$ example.cxx                             Example of use of matrix package
  2.  
  3. typedef double real;                 // all floating point double
  4.  
  5. #include "include.hxx"               // include standard files
  6.  
  7. #include "newmatap.hxx"              // need matrix applications
  8.  
  9. real t3(real);                       // round to 3 decimal places
  10.  
  11. // demonstration of matrix package on linear regression problem
  12.  
  13. main()
  14. {
  15.    // the data
  16.    // you may need to read this data using cin if you are using a
  17.    // compiler that doesn't understand aggregates
  18.    real y[]  = { 8.3, 5.5, 8.0, 8.5, 5.7, 4.4, 6.3, 7.9, 9.1 };
  19.    real x1[] = { 2.4, 1.8, 2.4, 3.0, 2.0, 1.2, 2.0, 2.7, 3.6 };
  20.    real x2[] = { 1.7, 0.9, 1.6, 1.9, 0.5, 0.6, 1.1, 1.0, 0.5 };
  21.    for (i=0; i<9;i++) cin >> y[i];
  22.    for (i=0; i<9;i++) cin >> x1[i];
  23.    for (i=0; i<9;i++) cin >> x2[i];
  24.    int nobs = 9;                           // number of observations
  25.    int npred = 2;                          // number of predictor values
  26.    int npred1 = npred+1;
  27.  
  28.    // we want to find the values of a,a1,a2 to give the best
  29.    // fit of y[i] with a0 + a1*x1[i] + a2*x2[i]
  30.  
  31.    // this example demonstrates three methods of calculation
  32.  
  33.    {
  34.       // traditional sum of squares and products method of calculation
  35.       // with subtraction of means
  36.  
  37.       // make matrix of predictor values
  38.       Matrix X(nobs,npred);
  39.  
  40.       // load x1 and x2 into X
  41.       //    [use << rather than = with submatrices and/or loading arrays]
  42.       X.Col(1) << x1;  X.Col(2) << x2;
  43.  
  44.       // vector of Y values
  45.       ColumnVector Y(nobs); Y << y;
  46.  
  47.       // make vector of 1s
  48.       ColumnVector Ones(nobs); Ones = 1.0;
  49.  
  50.       // calculate means (averages) of x1 and x2 [ .t() takes transpose]
  51.       RowVector M = Ones.t() * X / nobs;
  52.  
  53.       // and subtract means from x1 and x1
  54.       Matrix XC(nobs,npred);
  55.       XC = X - Ones * M;
  56.  
  57.       // do the same to Y [need (1,1) to convert 1x1 matrix to scalar]
  58.       ColumnVector YC(nobs);
  59.       real m = Matrix(Ones.t() * Y)(1,1) / nobs;  YC = Y - Ones * m;
  60.  
  61.       // form sum of squares and product matrix
  62.       //    [use << rather than = for copying Matrix into SymmetricMatrix]
  63.       SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  64.  
  65.       // calculate estimate
  66.       //    [bracket last two terms to force this multiplication first]
  67.       //    [ .i() means inverse, but inverse is not explicity calculated]
  68.       ColumnVector A = SSQ.i() * (XC.t() * YC);
  69.  
  70.       // calculate estimate of constant term
  71.       real a = m - Matrix(M * A)(1,1);
  72.  
  73.       // Get variances of estimates from diagonal elements of invoice of SSQ
  74.       //    [ we are taking inverse of SSQ; would have been better to use
  75.       //        CroutMatrix method - see documentation ]
  76.       Matrix ISSQ = SSQ.i(); DiagonalMatrix D; D << ISSQ;
  77.       ColumnVector V = D.CopyToColumn();
  78.       real v = 1.0/nobs + Matrix(M * ISSQ * M.t())(1,1);
  79.                         // for calc variance const
  80.  
  81.       // Calculate fitted values and residuals
  82.       ColumnVector Fitted = X * A + a;
  83.       ColumnVector Residual = Y - Fitted;
  84.       real ResVar = Residual.SumSquare() / (nobs-npred1);
  85.  
  86.       // print out answers
  87.       cout << "\nEstimates and their standard errors\n\n";
  88.       cout << a <<"\t"<< sqrt(v*ResVar) << "\n";
  89.       for (int i=1; i<=npred; i++)
  90.      cout << A(i) <<"\t"<< sqrt(V(i)*ResVar) << "\n";
  91.       cout << "\nObservations, fitted value, residual value\n";
  92.       for (i=1; i<=nobs; i++)
  93.      cout << X(i,1) <<"\t"<< X(i,2) <<"\t"<< Y(i) <<"\t"<<
  94.      t3(Fitted(i)) <<"\t"<< t3(Residual(i)) <<"\n";
  95.       cout << "\n\n";
  96.    }
  97.  
  98.    {
  99.       // traditional sum of squares and products method of calculation
  100.       // with subtraction of means - using Cholesky decomposition
  101.  
  102.       Matrix X(nobs,npred);
  103.       X.Col(1) << x1;  X.Col(2) << x2;
  104.       ColumnVector Y(nobs); Y << y;
  105.       ColumnVector Ones(nobs); Ones = 1.0;
  106.       RowVector M = Ones.t() * X / nobs;
  107.       Matrix XC(nobs,npred);
  108.       XC = X - Ones * M;
  109.       ColumnVector YC(nobs);
  110.       real m = Matrix(Ones.t() * Y)(1,1) / nobs;  YC = Y - Ones * m;
  111.       SymmetricMatrix SSQ; SSQ << XC.t() * XC;
  112.  
  113.       // Cholesky decomposition of SSQ
  114.       LowerTriangularMatrix L = Cholesky(SSQ);
  115.  
  116.       // calculate estimate
  117.       ColumnVector A = L.t().i() * (L.i() * (XC.t() * YC));
  118.  
  119.       // calculate estimate of constant term
  120.       real a = m - Matrix(M * A)(1,1);
  121.  
  122.       // Get variances of estimates from diagonal elements of invoice of SSQ
  123.       DiagonalMatrix D; D << L.t().i() * L.i();
  124.       ColumnVector V = D.CopyToColumn();
  125.       real v = 1.0/nobs + (L.i() * M.t()).SumSquare();
  126.  
  127.       // Calculate fitted values and residuals
  128.       ColumnVector Fitted = X * A + a;
  129.       ColumnVector Residual = Y - Fitted;
  130.       real ResVar = Residual.SumSquare() / (nobs-npred1);
  131.  
  132.       // print out answers
  133.       cout << "\nEstimates and their standard errors\n\n";
  134.       cout << a <<"\t"<< sqrt(v*ResVar) << "\n";
  135.       for (int i=1; i<=npred; i++)
  136.      cout << A(i) <<"\t"<< sqrt(V(i)*ResVar) << "\n";
  137.       cout << "\nObservations, fitted value, residual value\n";
  138.       for (i=1; i<=nobs; i++)
  139.      cout << X(i,1) <<"\t"<< X(i,2) <<"\t"<< Y(i) <<"\t"<<
  140.      t3(Fitted(i)) <<"\t"<< t3(Residual(i)) <<"\n";
  141.       cout << "\n\n";
  142.    }
  143.  
  144.    {
  145.       // Householder triangularisation method
  146.  
  147.       // load data - 1s into col 1 of matrix
  148.       Matrix X(nobs,npred1); ColumnVector Y(nobs);
  149.       X.Col(1) << 1.0;  X.Col(2) << x1;  X.Col(3) << x2;  Y << y;
  150.  
  151.       // do Householder triangularisation
  152.       // no need to deal with constant term separately
  153.       Matrix XT = X.t();             // Want data to be along rows
  154.       RowVector YT = Y.t();
  155.       LowerTriangularMatrix L; RowVector M;
  156.       HHDecompose(XT, L); HHDecompose(XT, YT, M); // YT now contains resids
  157.       ColumnVector A = L.t().i() * M.t();
  158.       ColumnVector Fitted = X * A;
  159.       real ResVar = YT.SumSquare() / (nobs-npred1);
  160.  
  161.       // get variances of estimates
  162.       L << L.i();
  163.       DiagonalMatrix D; D << L.t() * L;
  164.  
  165.       // print out answers
  166.       cout << "\nEstimates and their standard errors\n\n";
  167.       for (int i=1; i<=npred1; i++)
  168.      cout << A(i) <<"\t"<< sqrt(D(i)*ResVar) << "\n";
  169.       cout << "\nObservations, fitted value, residual value\n";
  170.       for (i=1; i<=nobs; i++)
  171.      cout << X(i,2) <<"\t"<< X(i,3) <<"\t"<< Y(i) <<"\t"<<
  172.      t3(Fitted(i)) <<"\t"<< t3(YT(i)) <<"\n";
  173.       cout << "\n\n";
  174.    }
  175. }
  176.  
  177. real t3(real r) { return int(r*1000) / 1000.0; }
  178.